/* Copyright 2007 dnAnalytics Project.
 *
 * Contributors to this file:
 * Patrick van der Velde
 * Marcus Cuda
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 */
using System;
using dnAnalytics.Math;
using dnAnalytics.Properties;

namespace dnAnalytics.LinearAlgebra
{
    public partial class DenseMatrix
    {
        /// <summary>
        /// Multiplies each element of this matrix with a scalar.
        /// </summary>
        /// <param name="scalar">The scalar to multiply with.</param>
        public override void Multiply(double scalar)
        {
            if (Precision.EqualsWithinDecimalPlaces(1.0, scalar, Constants.DecimalsForComparison))
            {
                return;
            }
            for (int i = 0; i < mData.Length; i++)
            {
                mData[i] *= scalar;
            }
        }

        /// <summary>
        /// Returns the transpose of this matrix.
        /// </summary>        
        /// <returns>The transpose of this matrix.</returns>
        public override Matrix Transpose()
        {
            DenseMatrix ret = new DenseMatrix(Columns, Rows);
            for (int j = 0; j < Columns; j++)
            {
                int index = j*Rows;
                for (int i = 0; i < Rows; i++)
                {
                    ret.mData[i*Columns + j] = mData[index + i];
                }
            }
            return ret;
        }

        /// <summary>
        /// Pointwise multiplies this matrix with another matrix and stores the result into the result matrix.
        /// </summary>
        /// <param name="other">The matrix to pointwise multiply with this one.</param>
        /// <param name="result">The matrix to store the result of the pointwise multiplication.</param>
        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null"/>.</exception>
        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null"/>.</exception>
        /// <exception cref="NotConformableException">If this matrix and <paramref name="other"/> are not the same size.</exception>
        /// <exception cref="NotConformableException">If this matrix and <paramref name="result"/> are not the same size.</exception>
        public override void PointwiseMultiply(Matrix other, Matrix result)
        {
            if (other == null)
            {
                throw new ArgumentNullException("other");
            }

            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (Columns != other.Columns || Rows != other.Rows)
            {
                throw new NotConformableException("result");
            }

            if (Columns != result.Columns || Rows != result.Rows)
            {
                throw new NotConformableException("result");
            }

            if (result is DenseMatrix && other is DenseMatrix)
            {
                DenseMatrix denseResult = (DenseMatrix) result;
                DenseMatrix denseOther = (DenseMatrix) other;
                for (int i = 0; i < mData.Length; i++)
                {
                    denseResult.mData[i] = mData[i]*denseOther.mData[i];
                }
            }
            else
            {
                base.PointwiseMultiply(other, result);
            }
        }

        /// <summary>
        /// Subtracts another matrix from this matrix.
        /// </summary>
        /// <param name="other">The matrix to subtract.</param>
        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null" />.</exception>
        /// <exception cref="NotConformableException">If the two matrices don't have the same dimensions.</exception>
        public override void Subtract(Matrix other)
        {
            if (other == null)
            {
                throw new ArgumentNullException("other");
            }

            if (other.Rows != Rows || other.Columns != Columns)
            {
                throw new NotConformableException("other", Resources.ParameterNotConformable);
            }

            if (other is DenseMatrix)
            {
                DenseMatrix otherDense = (DenseMatrix) other;

                for (int i = 0; i < mData.Length; i++)
                {
                    mData[i] -= otherDense.mData[i];
                }
            }
            else
            {
                base.Subtract(other);
            }
        }

        /// <summary>
        /// Adds another matrix to this matrix.
        /// </summary>
        /// <param name="other">The matrix to add to this matrix.</param>
        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null" />.</exception>
        /// <exception cref="NotConformableException">If the two matrices don't have the same dimensions.</exception>
        public override void Add(Matrix other)
        {
            if (other == null)
            {
                throw new ArgumentNullException("other");
            }

            if (other.Rows != Rows || other.Columns != Columns)
            {
                throw new NotConformableException("other", Resources.ParameterNotConformable);
            }
            if (other is DenseMatrix)
            {
                DenseMatrix otherDense = (DenseMatrix) other;

                for (int i = 0; i < mData.Length; i++)
                {
                    mData[i] += otherDense.mData[i];
                }
            }
            else
            {
                base.Add(other);
            }
        }


        /// <summary>
        /// An implementation of the BLAS GEMM routine: C:= alpha*op(A) * op(B) + beta*C.
        /// where op is to transpose the matrix or not, and C is this matrix.
        /// </summary>
        /// <param name="alpha">The alpha value.</param>
        /// <param name="beta">The beta value</param>
        /// <param name="transposeA">if set to <c>true</c>, transpose <paramref name="a"/>.</param>
        /// <param name="transposeB">if set to <c>true</c>, transpose <paramref name="b"/>.</param>
        /// <param name="a">the A matrix.</param>
        /// <param name="b">the B matrix.</param>
        public override void Gemm(double alpha, double beta, bool transposeA, bool transposeB, Matrix a, Matrix b)
        {
            if (a == null)
            {
                throw new ArgumentNullException("a");
            }
            if (b == null)
            {
                throw new ArgumentNullException("b");
            }

            if( !(a is DenseMatrix) || !(b is DenseMatrix))
            {
                base.Gemm(alpha, beta, transposeA, transposeB, a, b);
                return;
            }

            int aRows = a.Rows;
            int bRows = b.Rows;
            int aColumns = a.Columns;
            int bColumns = b.Columns;

            if (transposeA && transposeB)
            {
                if (aRows != bColumns)
                {
                    throw new NotConformableException();
                }
                if (Rows != aColumns && Columns != bRows)
                {
                    throw new NotConformableException();
                }
            }
            else if (transposeA)
            {
                if (aRows != bRows)
                {
                    throw new NotConformableException();
                }
                if (Rows != aColumns && Columns != bColumns)
                {
                    throw new NotConformableException();
                }
            }
            else if (transposeB)
            {
                if (aColumns != bColumns)
                {
                    throw new NotConformableException();
                }
                if (Rows != aRows && Columns != bRows)
                {
                    throw new NotConformableException();
                }
            }
            else
            {
                if (aColumns != bRows)
                {
                    throw new NotConformableException();
                }
                if (Rows != aRows && Columns != bColumns)
                {
                    throw new NotConformableException();
                }
            }

            if (Precision.EqualsWithinDecimalPlaces(0.0, alpha, Constants.DecimalsForComparison) && Precision.EqualsWithinDecimalPlaces(0.0, beta, Constants.DecimalsForComparison))
            {
                Clear();
                return;
            }

            double[] adata;
            if (ReferenceEquals(a, this) || ReferenceEquals(a, b))
            {
                adata = new DenseMatrix(a).mData;
            }
            else
            {
                adata = ((DenseMatrix) a).mData;
            }

            double[] bdata;
            if (ReferenceEquals(a, this) || ReferenceEquals(a, b))
            {
                bdata = new DenseMatrix(b).mData;
            }
            else
            {
                bdata = ((DenseMatrix) b).mData;
            }
            double[] cdata = mData;

            if (Precision.EqualsWithinDecimalPlaces(1.0, alpha, Constants.DecimalsForComparison))
            {
                if (Precision.EqualsWithinDecimalPlaces(0.0, beta, Constants.DecimalsForComparison))
                {
                    if (transposeA && transposeB)
                    {
                        for (int j = 0; j != aColumns; j++)
                        {
                            int jIndex = j*Rows;
                            for (int i = 0; i != bRows; i++)
                            {
                                int iIndex = i*aRows;
                                double s = 0;
                                for (int l = 0; l != bColumns; l++)
                                {
                                    s += adata[iIndex + l]*bdata[l*bRows + j];
                                }
                                cdata[jIndex + i] = s;
                            }
                        }
                    }
                    else if (transposeA)
                    {
                        for (int j = 0; j != bColumns; j++)
                        {
                            int jcIndex = j*Rows;
                            int jbIndex = j*bRows;
                            for (int i = 0; i != aColumns; i++)
                            {
                                int iIndex = i*aRows;
                                double s = 0;
                                for (int l = 0; l != aRows; l++)
                                {
                                    s += adata[iIndex + l]*bdata[jbIndex + l];
                                }
                                cdata[jcIndex + i] = s;
                            }
                        }
                    }
                    else if (transposeB)
                    {
                        for (int j = 0; j != bRows; j++)
                        {
                            int jIndex = j*Rows;
                            for (int i = 0; i != aRows; i++)
                            {
                                double s = 0;
                                for (int l = 0; l != aColumns; l++)
                                {
                                    s += adata[l*aRows + i]*bdata[l*bRows + j];
                                }
                                cdata[jIndex + i] = s;
                            }
                        }
                    }
                    else
                    {
                        for (int j = 0; j != bColumns; j++)
                        {
                            int jcIndex = j*Rows;
                            int jbIndex = j*bRows;
                            for (int i = 0; i != aRows; i++)
                            {
                                double s = 0;
                                for (int l = 0; l != aColumns; l++)
                                {
                                    s += adata[l*aRows + i]*bdata[jbIndex + l];
                                }
                                cdata[jcIndex + i] = s;
                            }
                        }
                    }
                }
                else
                {
                    if (transposeA && transposeB)
                    {
                        for (int j = 0; j != aColumns; j++)
                        {
                            int jIndex = j*Rows;
                            for (int i = 0; i != bRows; i++)
                            {
                                int iIndex = i*aRows;
                                double s = 0;
                                for (int l = 0; l != bColumns; l++)
                                {
                                    s += adata[iIndex + l]*bdata[l*bRows + j];
                                }
                                cdata[jIndex + i] = cdata[jIndex + i]*beta + s;
                            }
                        }
                    }
                    else if (transposeA)
                    {
                        for (int j = 0; j != bColumns; j++)
                        {
                            int jcIndex = j*Rows;
                            int jbIndex = j*bRows;
                            for (int i = 0; i != aColumns; i++)
                            {
                                int iIndex = i*aRows;
                                double s = 0;
                                for (int l = 0; l != aRows; l++)
                                {
                                    s += adata[iIndex + l]*bdata[jbIndex + l];
                                }
                                cdata[jcIndex + i] = s + cdata[jcIndex + i]*beta;
                            }
                        }
                    }
                    else if (transposeB)
                    {
                        for (int j = 0; j != bRows; j++)
                        {
                            int jIndex = j*Rows;
                            for (int i = 0; i != aRows; i++)
                            {
                                double s = 0;
                                for (int l = 0; l != aColumns; l++)
                                {
                                    s += adata[l*aRows + i]*bdata[l*bRows + j];
                                }
                                cdata[jIndex + i] = s + cdata[jIndex + i]*beta;
                            }
                        }
                    }
                    else
                    {
                        for (int j = 0; j != bColumns; j++)
                        {
                            int jcIndex = j*Rows;
                            int jbIndex = j*bRows;
                            for (int i = 0; i != aRows; i++)
                            {
                                double s = 0;
                                for (int l = 0; l != aColumns; l++)
                                {
                                    s += adata[l*aRows + i]*bdata[jbIndex + l];
                                }
                                cdata[jcIndex + i] = s + cdata[jcIndex + i]*beta;
                            }
                        }
                    }
                }
            }
            else
            {
                if (transposeA && transposeB)
                {
                    for (int j = 0; j != aColumns; j++)
                    {
                        int jIndex = j*Rows;
                        for (int i = 0; i != bRows; i++)
                        {
                            int iIndex = i*aRows;
                            double s = 0;
                            for (int l = 0; l != bColumns; l++)
                            {
                                s += adata[iIndex + l]*bdata[l*bRows + j];
                            }
                            cdata[jIndex + i] = cdata[jIndex + i]*beta + alpha*s;
                        }
                    }
                }
                else if (transposeA)
                {
                    for (int j = 0; j != bColumns; j++)
                    {
                        int jcIndex = j*Rows;
                        int jbIndex = j*bRows;
                        for (int i = 0; i != aColumns; i++)
                        {
                            int iIndex = i*aRows;
                            double s = 0;
                            for (int l = 0; l != aRows; l++)
                            {
                                s += adata[iIndex + l]*bdata[jbIndex + l];
                            }
                            cdata[jcIndex + i] = alpha*s + cdata[jcIndex + i]*beta;
                        }
                    }
                }
                else if (transposeB)
                {
                    for (int j = 0; j != bRows; j++)
                    {
                        int jIndex = j*Rows;
                        for (int i = 0; i != aRows; i++)
                        {
                            double s = 0;
                            for (int l = 0; l != aColumns; l++)
                            {
                                s += adata[l*aRows + i]*bdata[l*bRows + j];
                            }
                            cdata[jIndex + i] = alpha*s + cdata[jIndex + i]*beta;
                        }
                    }
                }
                else
                {
                    for (int j = 0; j != bColumns; j++)
                    {
                        int jcIndex = j*Rows;
                        int jbIndex = j*bRows;
                        for (int i = 0; i != aRows; i++)
                        {
                            double s = 0;
                            for (int l = 0; l != aColumns; l++)
                            {
                                s += adata[l*aRows + i]*bdata[jbIndex + l];
                            }
                            cdata[jcIndex + i] = alpha*s + cdata[jcIndex + i]*beta;
                        }
                    }
                }
            }
        }
    }
}